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Abstract: The channeling of the ion recoihng after a colhsion with a WIMP changes the 
ionization signal in direct detection experiments, producing a larger signal than otherwise 
expected. We give estimates of the fraction of channeled recoiling ions in Nal (Tl) crystals 
using analytic models produced since the 1960's and 70's to describe channeling and blocking 
effects. We find that the channeling fraction of recoiling lattice nuclei is smaller than that of 
ions that are injected into the crystal and that it is strongly temperature dependent. 
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1. Introduction 

Ion channeling in crystals has received a large amount of attention in the interpretation of 
experiments designed to search for Weakly Interacting Massive dark matter Particles (dark 
matter WIMPs) through their scattering in a low-background detector. Channeling would 
occur when the nucleus that recoils after being hit by a dark matter particle moves off in a 
direction close to a symmetry axis or symmetry plane of the crystal. Channeled ions loose 
their energy predominantly to electrons, while non-channeled ions transfer their energy to 
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lattice nuclei. In scintillators like NaT (Tl), which are sensitive to the electronic energy losses, 
channeling increases the fraction of recoil energy that is observed as scintillation light. The 
presence of channeling in NaT (Tl), for example, affects the mass and cross section of the 
dark matter particles inferred from the observation of an annual modulation by the DAMA 
collaboration. As a consequence, the comparison between the DAMA results and the null 
results of other experiments is affected too . 

In addition, ion channeling in crystals could give rise to a daily modulation due to the 
preferred direction of the dark matter flux arriving on the Earth (this is the direction of the 
motion of the Solar System with respect to the Galaxy and its dark halo). Earth's daily 
rotation naturally changes the direction of the "WIMP wind" with respect to the crystal 
axes, thus the amount of recoiling ions that are channeled vs non-channeled, and therefore 
the amount of energy visible via scintillation. Avignone, Creswick, and Nussinov Q pointed 
this out for Nal (Tl) crystals, although the available estimates of the strength of these daily 
modulations are somewhat simplistic. 

Given the importance of channeling in the interpretation of direct detection experiments, 
and the need to refine the current calculations of a possible daily modulation due to chan- 
neling, it is worthwhile to take a deeper look at channeling in the context of dark matter 
detection. We set out to improve the evaluation of the channeling daily modulation, extend 
it from Nal (Tl) to Si, Ge and other crystals, and in the process understand channeling in 
dark matter experiments through analytical means. In this paper, the first one of a series, 
we present our analytic calculations of the channeling fraction in Nal (Tl). In the subsequent 
papers, we will present our results for Si, Ge, and other crystals, and our conclusions on the 
daily modulation due to channeling. 

2. Ion channeling and blocking 

Channeling and blocking effects in crystals refer to the orientation dependence of charged 
ion penetration in crystals. In the "channeling effect" ions incident upon a crystal along 
symmetry axes and planes suffer a series of small-angle scatterings that maintain them in 
the open "channels" between the rows or planes of lattice atoms and thus penetrate much 
further into the crystal than in other directions. Channeled incoming ions do not get close 
to lattice sites, where they would be deflected at large angles. The "blocking effect" consists 
in a reduction of the flux of ions originating in lattice sites along symmetry axes and planes, 
creating what is called a "blocking dip" in the flux of ions exiting from a thin enough crystal 
as a function of the exit angle with respect to a particular symmetry axis or plane. Directional 
effects in ion penetration in crystals was first observed in 1960 |^] in the sputtering ratio for 
ions bombarding a single crystal and its explanation in terms of channeling was first done 
in 1962 Q, although the effect had been predicted to exist in 1912 Strongly anisotropic 
effects for positive particle trajectories originating at lattice sites were discovered in 1965, with 
particles emitted from radioactive atoms and wide-angle scattering of positive ions in several 
experiments. Immediately the relation between these blocking effects and the channeling 
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effect was explained by Lindhard Q in 1965. In the 1960's and 70's the experimental and 
theoretical work on channeling proceeded at a very fast pace (see for example the review by 
D. Gemmell 0] and references therein). 

Channeling and blocking effects in crystals are used in crystallography, in the study of 
lattice disorder, ion implantation, and the location of dopant and impurity atoms in crys- 
tals, in studies of surfaces, interfaces and epitaxial layers, in measurements of short nuclear 
lifetimes, in the production of polarized beams etc (see for example @, |^). Channeling and 
blocking effects are related because the non-channeled incident ions are those which suffer 
a close-encounter process with an atomic nucleus in the crystal, namely those which pass 
sufficiently close to a lattice nucleus to be deflected at a large angle. After a close-encounter 
collision the deflected ion acts as if it was "emitted" from a lattice site. Channeling is many 
times observed as a lack of large angle deflections for ions incident at a small angle ■0 with 
respect to a particular symmetry axis or plane. This forms a "channeling dip" in the out- 
going flux as a function of the incident beam angle ip. As first pointed out by Lindhard 
when no slowing-down processes are involved, the "channeling" and "blocking" dips should 
be identical, when compared for the same particles, energies, crystals and crystal directions. 

Ion channeling in Nal (Tl) was first observed in 1973 by Altman, Dietrich, Murray and 
Rock [|l^]. They observed that the scintillation output of a monochromatic 10 MeV ^^O beam 
through an Nal (Tl) scintillator shows two peaks: one at low energy due to non-channeled 
ions, and one at high energy due to channeled ions. The channeled ions produce more 
scintillation light because they lose most of their energy via electronic stopping rather than 
nuclear stopping. 

This may be an important effect in direct dark matter detection experiments in which a 
scintillation signal due to the recoil of ions as a result of WIMP collisions is searched for. The 
potential importance of the channeling effect for direct dark matter detection was first pointed 
out for Nal (Tl) by Drobyshevski [11 1 and by the DAMA collaboration [12|. When Na or I 
ions recoiling after a collision with a dark matter WIMP move along crystal axes and planes, 
their quenching factor is approximately Q = 1 instead of QNa = 0.3 and Qi = 0.09, since 
they give their energy to electrons. The DAMA collaboration |12] found that the fraction of 
channeled recoils is large for low recoiling energies in the keV range, and that this effect shifts 
the regions in cross section versus mass of acceptable WIMP models in agreement with the 
DAMA data towards lower WIMP masses. 

Most of the applications of channeling and blocking are at energies of MeV and higher, 
however some use much lower energies, up to the keV range. In particular, avoiding chan- 
neling is essential in the manufacturing of semiconductor devices, since ion implantation at a 
controlled depth is the primary technique. Boron, arsenic and phosphorus ions are implanted 
in silicon, for example, to produce integrated circuits, at energies from lOO's of eV to several 
MeV (see for example [p^]). 

In this paper we compute the channeling fraction of recoiling ions in Nal (Tl) crystals as 
function of the recoil energy and the temperature. We do the same for Ge and Si crystals in 
a companion paper. 
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3. Models of Channeling 



3.1 Continuum models 

There are different approaclies to calculate tlie deflections of ions traveling in a crystal. In 
the "binary collision model" the ion path is computed by a computer program (see Ref. [T^ 
for one of the first ones) in terms of a succession of individual interactions, each with one 
of the atoms in the crystal. Crystal imperfections and lattice vibrations are thus easily and 
correctly taken into account. In "continuum models", reasonable approximations are made 
which allow to replace the discrete series of binary collisions with atoms by a continuous 
interaction between a projectile and uniformly charged strings or planes. These models allow 
to replace the numerical calculations by an analytic description of channeling, and provide 
good quantitative predictions of the behavior of projectiles in the crystal in terms of simple 
physical quantities. This is the approach we use here. This analytical description was initially 
developed mostly by J. Lindhard and collaborators for ions of energy MeV and higher, 
and its use was later extended to lower energies, i.e. hundreds of eV and above, mostly to 
apply it to ion implantation in Si. This approach must be complemented by determination 
of parameters through data fitting or simulations. Moreover, lattice vibrations are more 
difficult to include in continuum models. Since we use a continuum model, our results should 
in last instance be checked by using some of the many sophisticated simulation programs that 
implement the binary collision approach or mixed approaches (e.g. |15, |l^, |l7|]). 



Although the analytical description works better at higher energies (where it has been 
very well tested experimentally), at low and intermediate energies the critical angles for 
channeling predicted by analytic models have also been found to be in good agreement with 
experimental results. For the low energy range we found most useful the work of G. Hobler, 



who in 1995 and 1996 |18| perfected and checked experimentally previous continuum model 
predictions |19| for axial and planar channeling at energies in the keV to a few 100 keV range, 
to avoid channeling in the implantation of B, P and As atoms in Si crystals. Measurements of 
axial critical angles obtained in the late 1960's for light (H"*", D"*" and He"*") and intermediate 
mass (B and Ar) ions propagating in various crystals (gold, tungsten, silicon) with energies 
between 1 and 100 keV were found to be in good agreement with the predictions of Lindhard 



models 21, 22]. In 1999 K. M. Lui |2^ and collaborators compared experimental results 



of 5 keV Ne"*" ions on platinum and predictions of the trajectory simulation code SARIC |16], 
based on the binary collision model, with the predictions of Lindhard's analytical model and 
the observed angular half-width of the blocking dips for axial channels were found to be in 
good qualitative agreement with Lindhard's critical angle (both are similar as can be seen 
in Table 1 of Ref In 2002, S. M. Hogg et al. [||] studied channeled implantation of 

80 keV Er ions into Si and concluded that the axial measured critical angle was in excellent 
agreement with both computer simulations (made with the MDRANGE program ||T^) and 
experimental results. In 2005 Lindhard's critical angle prediction was used to understand 
qualitative features of computational results of the SARIC program for 4 keV Ne^ ions 



impinging on a Pt surface [25|. 
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Our calculation is based on the classical analytic models developed in the 1960's and 
70's, in particular by Lindhard 2C, 28, |3^, 31, The fact that the de Broglie 



wavelengths of ions in the keV-MeV range are of the order of ~ 0.01 pm (and smaller at 
higher energies), which is much less than the lattice constant of a crystal (~ 10 pm), justifies 
using a classical treatment. We use the continuum string and plane model, in which the 
screened Thomas-Fermi potential is averaged over a direction parallel to a row or a plane. 
This averaged potential U is considered to be uniformly smeared along the row or plane of 
atoms, which is a good approximation if the propagating ion interacts with many lattice 
atoms in the row or plane by a correlated series of many consecutive glancing collisions with 
lattice atoms. We are going to consider just one row, which simplifies the calculations and is 
correct except at the lowest energies we consider, as we explain below. 

There are several good analytic approximations of the screened potential. Here we use 
Lindhard's expression, because it is the simplest and allows to find analytical expressions for 
the quantities we need. The transverse averaged continuum potential of a string as a function 
of the transverse distance r to the string, relevant for axial channeling, was approximated by 
Lindhard |^] as 

U{r) =E^Pl -In + (3.1) 

where C is a constant, which was found experimentally to be C ~ \/3, and 

^1 = (3.2) 

Zi , Z2 are the atomic numbers of the recoiling and lattice nuclei respectively, d is the spacing 
between atoms in the row, a is the Thomas-Fermi screening distance, a = 0.4685A(Z]^^^ -|- 
Z2^'^)~'^/^ |14, |7| and E = Mv'^/2 is the kinetic energy of the propagating ion. In our case. 



E is the recoil energy imparted to the ion after a collision with a WIMP, 

E=--^, (3.3) 
2M' ^ ^ 

where q is the recoil momentum. The string of crystal atoms is at r = 0. 

The transverse averaged continuum potential of a plane of atoms, relevant for planar 
channeling, given by Lindhard |^ as a function of the distance x perpendicular to the plane 
is 



J a 



(3.4) 



where ipa is 

( 

and n = Ndpch is the average number of atoms per unit area, where is the atomic density 
and dpch is the width of the planar channel, i.e. the interplanar spacing (thus the average 



i^a = ( ^""^f " ) (3.5) 
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Figure 1: Continuum axial (black) and planar (green/gray) potentials for (a) Na and (b) I ions, 
propagating in the <100> axial and {100} planar channels of an Nal crystal. The screening radii 
shown as vertical lines are ONa = 0.00878 nm and a/ — 0.0115 nm (see App. A). 



distance of atoms within a plane is dp = 1/ y/Ndpch)- The plane is at a; = 0. Examples of 
axial and planar continuum potentials are shown in Fig. |^. 

The continuum model does not imply that the potential energy of an ion moving near 
an atomic row is well approximated by the continuum potential U. The actual potential 
consists of sharp peaks near the atoms and deep valleys in between. The continuum model 
says that the net deflection due to the succession of impulses from the peaks is identical to 
the deflection due to a force —U'. This is only so if the ion never approaches any individual 
atom so closely that it suffers a large-angle collision. Lindhard proved that for a string of 
atoms this is so only if 

U"{r) < ^E, (3.6) 

where the double prime denotes the second derivative with respect to r. Replacing the 
inequality in Eq. 3.6 by an equality defines an energy dependent critical distance Tc such 
that r > Tc for the continuum model to be valid. Morgan and Van Vliet ||2^ also derived a 
condition for axial channels, similar to Eq. (but with the factor 8 replaced by 16). 

The condition in Eq. |3.6| for the validity of the continuum model on the axial effective 
potential is equivalent (as shown initially by Lindhard and proven below) to insuring that the 
minimum distance of approach to the string remains larger than ~ d for large E (MeV 
and above) and ~ d -02 (with ip2 given in Eq. 3.1^ ) for small E (below lOO's of keV). Thus, 
the smaller the atomic interdistance d and the larger the ion velocity (i.e. the smaller ipi or 
■02) the more accurate the continuum model Q. 

The appearance of the angle '01 in this condition can be easily understood by considering 
the "Coulomb shadow" formed by individual atoms behind the direction of arrival of a parallel 
beam of positive ions. For an unscreened Coulomb potential and small angle deflections, the 
trajectories of the projectiles give rise to a shadow cone in which the projectiles do not enter. 
A distance z behind the deflecting nucleus in the direction of arrival of the projectiles, the 
shadow cone radius is D{z) = \j2zd 11)1 (here d enters through the definition of 0i). As the 



- 6 - 



incident angle of the incoming ions with respect to the row of atoms decreases, there is a 
critical value of the incident angle at which the edge of the shadow of an atom passes through 
the adjacent atom. This critical angle is approximately D{d)/d = \pl tpi- For angles of 
incidence larger than ~ ipi , the shadow cones of the atoms in a row are independent of each 
other. But for incident angles smaller than ~ ^i, the shadow cones interfere with each other, 
so that the atoms in a row are effectively shadowed and not exposed to the projectiles ^. 
In this case the incident ions do not approach the shadowed atoms closer than a distance 
D{d) ~ d ipi, as mentioned above. For a more realistic screened Coulomb potential (as 
considered in this paper) the shadow cone radius is smaller than for the unscreened potential 
and the difference between the two becomes smaller for higher ion energies. 

The breakdown of the continuum theory for planar channeling is more involved than 
for axial channeling because the atoms in the plane contributing to the scattering of the 
propagating ion are usually displaced laterally within the plane. Thus the moving ion does 
not encounter atoms at a fixed separation or at fixed impact parameter as is the case for a 
row. Morgan and Van Vliet reduced the problem of scattering from a plane of atoms to 
the scattering from an equivalent row of atoms contained in a strip centered on the projection 



of the ion path on the plane of atoms. They then applied Eq. 3.6 to the fictitious string 



defined in this way as the condition for planar channeling (more about this below). 
3.2 The transverse energy 

Lindhard proved that for channeled particles the longitudinal component vcos<j), i.e. the 
component along the direction of the row or plane of the velocity, may be treated as constant 
(if energy loss processes are neglected). Then, in the continuum model, the trajectory of 
the ions can be completely described in terms of the transverse direction, perpendicular to 
the row or plane considered. For small angle between the ion's trajectory and the atomic 
row (or plane) in the direction perpendicular to the row (or plane), the so called "transverse 
energy" 

E± = E sin^ (p + U Ecp"^ + U (3.7) 



is conserved. In Eq. 3.7 relativistic corrections are neglected. 

In each binary collision of the ion with the closest atom, E± changes abruptly, because 
the angle (j) changes in a very short time (i.e. while the potential is practically constant). 
Then, between two collisions, the change is compensated because the potential component of 
E± changes continuously as the ion propagates while the angle cp is constant . A good way 
to test to what extent this compensation takes place is to calculate the value of E± far away 
from the collision sites, namely half-way between successive collision sites. 

The condition in Eq. |3.6| was derived by Lindhard (in appendix A of Ref. [^) for axial 
channels by defining E± at the planes half-way between string of atoms and asking for E± at 
contiguous half-way planes to be conserved to first order (this is the so called "half-way plane" 



model). Morgan and Van Vliet |27] derived a condition for axial channeling very similar to 



Eq. |3.6| by calculating the difference between the scattering angle due to a binary collision 
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and the deflection angle in the continuum potential when the ion travels the distance between 
two contiguous halfway planes (the half-way planes considered by Lindhard). 

Let be the initial position at which the WIMP nucleus collision occurs, i.e. if > the 
recoiling nucleus was displaced with respect to its position of equilibrium in the string when 
it collided with a WIMP. We call (j) the angle of the initial recoil momentum with respect to 
the row of atoms, and E the initial recoil energy of the propagating ion. Given these initial 
parameters, the issue of where to define Ex_ arises. Namely, we define 

E^ = EW ^ + U{r*), (3.8) 

but there are different possible choices for r*, the position at which to measure the potential 
U . In the "half-way plane model" used by Lindhard, U is measured after the recoiling ion 
propagates a distance d/2 along the string, when it is at a distance 

r* = rJip = n + {d/2) tan (j) (3.9) 

perpendicular to the string at the halfway-plane. All angles we are dealing with are small 
enough that svncj) ~ ~ cf). This choice was shown to work better in some respects [2S] 

(such as the blocking angular distribution in axial channels) than the "continuum approxi- 
mation." In the latter, the transverse energy E± is considered to be conserved all along the 
string, not only at the halfway-planes, in which case r* is chosen to be just 

r* = r*cA = n. (3.10) 

The two choices r* = r^p and r* = coincide only if dtan(/)/2rj ^ 1, a condition that 
at energies below lOO's of keV is in general not fulfilled, in which case the "continuum ap- 
proximation" is not a good approximation. In fact, assuming the "continuum approxima- 
tion", the angle at the first halfway plane must have a different value, 0' say, such that 
sin^ 0' — sin^ (/) = [[/(r^) — [/(rj + ^dtan (/))]/£^. Fig. |l] shows that the potential C/ at a distance 
a or larger is in the 1-10 keV range. Thus the difference between both definitions of the 
transverse energy is very small at large enough values of the energy ^ 10 keV, but for 
lower values of E^ in the keV to the lO's of keV range, the definitions r* = rjjp and r* = 
give different results unless (j) is small enough. In Ref. |^9| the predictions of both models 
for blocking of 400 keV protons in W and 7 MeV protons in Si were compared with the 
predictions of the binary-collision model. The "half-way plane model" results were found to 
be in agreement with those of the binary-collision model, even when those of the "continuum 
approximation" were not. 

In all these cases, even when considering blocking, the propagating ion was always differ- 
ent than a lattice ion. In our case, the recoiling ion leaves an empty lattice site, thus it moves 
away from an empty lattice site in the potential generated by its neighboring lattice atoms. 
So the potential the recoiling ion moves through at the moment of collision is very small, and 
the recoiling ion conserves its momentum and direction of motion until it gets very near the 
nearest neighbor, a distance d away along the string. At this moment, it is at a distance 

r* = r*gj. = r j + dtan^j (3-11) 
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from its nearest neighbor. Therefore, we will make the approximation of defining the potential 



entering into Eq. 3.8 at r* = r 



3.3 Minimum distance of approach and critical channeling angle 

The conservation of the transverse energy provides a definition of the minimum distance of 
approach to the string, rmm (or to the plane of atoms Xmin), at which the trajectory of the 
ion makes a zero angle with the string (or plane), and also of the angle ifj at which the ion 
exits from the string (or plane), i.e. far away from it where [/ ~ 0. In reality the furthest 
position from a string or plane of atoms is the middle of the channel, whose width we call 
dach for an axial channel or d^ch for a planar channel, thus 

= U{r^i^) = E^^ + C/(dach/2). (3.12) 

We define the axial channel width dach in terms of the interatomic distance d as dach = 1/ V Nd, 
where is the atomic density. 

For axial channeling Lindhard equates the condition for channeling with the condition in 
Eq. |3.6| for the validity of the continuum model. For Lindhard's axial potential, this condition 
reads 

^ Eid'^ l + 3(^f|if)2 ^ ^ 



r: 



2 a + (1^)2 



mm 



where Ei = Eipf (and ipi was defined in Eq. |3.2| ). Since the right-hand side of this inequality is 
a monotonically decreasing function of r^in, one just needs to solve the equation obtained by 
replacing the inequality with an equality. Solving a cubic equation, the condition in Eq. |3.13 



can be inverted to find rc{E), the minimum value of rmin- We find the following decreasing 
function of E 



r,{E)=Ca^ 
where to simplify the expression we defined 



, /I (l-3z/2)\ 

Vl + cos - arccos — —j^r - 1 

3 (1 + ^)3/2 J 
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•2 



(3.14) 



(3.15) 



SEC^a^' 

This expression gives the smallest possible minimum distance of approach of the propa- 
gating ion with the row for a given energy E, i.e. r^i^ > r(.{E) and, since the potential U (r) 
decreases monotonically with increasing r, 

?7(rmm) < U{r,{E)). (3.16) 



Using Eq. this can be further translated into an upper bound on E±_ and on -0, the angle 



the ion makes with the string far away from it. 
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ipc{E) is the maximum angle the ion can make with the string far away from it (i.e. in the 
middle of the channel) if the ion is channeled. When C/((iach/2) can be neglected, i.e. when 
rc{E) < (iach/2, the limiting values of ipdE) (as already proven by Lindhard p) are i^c{E) ~ 



ipi (see Eq. 3.2) for large E (z 1, typically close to MeV and larger) and ipc{E) — ip2 at 



low E {z ^ 1, typically smaller than a few 100 keV), where 



*'-f^- (3..) 

One can easily see that the critical distance Tc becomes Vc — d ipi/2\/2 for large E and 
rc — d ip2 for small E. 

The critical distance rc{E) increases as E decreases. At low enough E, rc{E) becomes 
close to (iach/2, and the critical angle ipc{E), the maximum angle for channeling in the middle 
of the channel, goes to zero. This means that there is a minimum energy below which 
channeling cannot happen, even for ions moving initially in the middle of the channel. This is 
a reflection of the fact that the range of the interaction between ion and lattice atoms increases 
with decreasing energy and at some point there is no position in the crystal where the ion 
would not be deflected at large angles. The existence of a minimum energy for channeling 
was found by Rozhkov and Dyuldya |^4| in 1984 and later by Hobler [18| in 1996. It is clear 



that to compute rc{E) when it is not small with respect to (iach/2, and thus to compute the 
actual minimum energy for channeling, we would need to consider the effect of more than one 
row or plane (as done in Refs. [Q] and jl^]), thus our results are approximate in this case. 
For planar channeling we will follow the definition of fictitious row in Morgan and Van 



Vliet 1 27, 18 1 . They reduced the problem of scattering from a plane of atoms to the scattering 
of an equivalent row of atoms contained in a strip of width 2R (R is defined below) centered 
on the projection of the ion path on the plane of atoms, and took the average area per atom 
in the plane to be 2R times the characteristic distance d between atoms along this fictitious 
row, 

d=l/{Ndp,^2R). (3.19) 
Once the width 2R of the fictitious row is specified, one uses the channeling condition for the 



continuum string model, Eq. 3.6, with an average atomic composition of the plane. For R, 



Morgan and Van Vliet used the impact parameter in an ion-atom collision corresponding to a 
defiection of the order of the break-through angle y/Up{0)/E. This is the defiection necessary 
for an ion of energy E approaching the plane from far away (so that the initial potential can 
be neglected) to overcome the potential barrier at the center of the plane at x = (namely 
so that Ej_ = Up{0)). Using the Moliere approximation for the screened potential (which we 
do not use), Morgan and Van Vliet found for d 

-1 



d 



-MV 



A aiVdpchln B ZiZ2e^ / a\ EUp{d) 



(3.20) 



with coefficients A = 1.2 and B = 4. However, Morgan and Van Vliet [27| found discrepancies 



with this theoretical formula in simulations of binary collisions of 20 keV protons in a copper 
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crystal and adjusted the coefficients to yl = 3.6 and B = 2.5. Hobler ||18[] used both sets of 
coefficients and compared them with simulations and data of B and P in Si for energies of 
about 1 keV and above. Hobler concluded that the original theoretical formula was better in 
his case. In any case, Hobler proposed yet another empirical relation to define d. 

In the absence of simulations for Nal, we are going to use an upper bound on i?, given 
by the average interdistance of atoms in the plane, 2R < dp = 1/ -s/ Ndpch, so that replacing 
the maximum value of R in Eq. 3.19| we find that the minimum value of d is the average 



interdistance of atoms in the plane, dp 

drain — dp. (3.21) 



Thus, for planar channeling, we use the condition in Eq. 3.6 for a fictitious string, re- 
placing the distance d by the distance dp and replacing the composition of the string for the 



average composition of the plane. Let us call rc{E) the critical distance obtained from Eq. p. 14 



for this fictitious string, then the minimum distance of approach for planar channeling is 

Xc{E) = r^{E). (3.22) 

The use of d = dp yields a lower bound on Xc, as shown in Fig. |2| (and thus an upper bound on 
the fraction of channeled recoils as explained later) . Fig. ^ shows the planar critical distances 
of approach Xc (Eq. |3.22| ) using the theoretical (with coefficients A = 1.2 and B = 4) and 
adjusted (with coefficients A = 3.6 and B = 2.5) Morgan-Van Vliet expressions for d in 
Eq. p. 20 . Our choice of Xc, with d in Eq. [3.21 is also plotted in Fig. ^. We can see that it is 



lower than the others. 

Writing equations equivalent to Eq. 3.12 and |3.1(: for planar channels, namely 



El. = U{xmrn) = E{ipPf + C/p(dpch/2), (3.23) 

and 

C/p(Xmm) < Up{x,{E)), (3.24) 
we obtain an equation similar to Eq. |3.17| but for the maximum planar channeling angle. 



rAE) = ^MffM^MW?). (3.25) 

For very small energies, for which Xc{E) > dpch/2 no channeling is possible (the maximum 
distance to any plane cannot be larger than half the width of the channel separating them) . 
When Xc{E) approaches the middle of the channel the effect of other planes should be con- 
sidered, so our approximation of using the potential of only one plane is not correct in this 
regime. 

The static lattice critical distances are presented in left panel of Figs. |^ and ^ for the 
100 and 111 axial and planar channels. 

There is an alternative way of treating planar channels presented by Matyukhin | |35[ | in 
2008, but we have doubts about the validity of this method, for which we have not found any 
comparison with either simulations or data. For completeness we present it in Appendix C. 
It predicts larger channeling fractions. 
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Figure 2: Comparison of planar critical distances of approach as given by the theoretical (solid line) 
and adjusted (dashed line) Morgan- Van Vliet (label MV, green/gray) expressions for d and our choice 
(black) of Xc, which we take as a lower bound, as function of the energy for Na ions travelling in {100} 
planar channels. Also shown is the radius of the channel (ipch/2. 



3.4 Temperature dependent critical distances and angles 

So far we have been considering static strings and planes, but the atoms in a crystal are actu- 
ally vibrating. We use here the Debye model to take into account the zero point energy and 
thermal vibrations of the atoms in a crystal. The one dimensional rms vibration amplitude 
ui of the atoms in a crystal in this model is 0, ^] 



ui{T) = 12.1 A 



*<®''^' ' ^1 (Me)- 



1/2 



(3.26) 



e/T 4^ 

where the 1/4 term accounts for the zero point energy, M is the atomic mass in amu, Q and 
T are the Debye temperature and the temperature of the crystal in K, respectively, and $(x) 
is the Debye function, 

X Jo e* - 1 

Eq. 3.26| was derived for monoatomic cubic crystals for which M is clearly specified. In the 
case of crystals composed of more than one kind of atom, experiments have shown that the 
difference of vibration amplitudes of both types is very small for T > Q [^ 1, even when the 
difference of atomic weights of the various kinds of atoms is large, as is the case for Nal. Using 
as M the average mass 

M = (MNa + Mi)/2 (3.28) 

produces an error of less than 10% in the actual vibration amplitudes at T > p6|| . For Nal, 
we take the Debye temperature to be = 165 K [0, ^ (although it changes with T between 
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169 K at a few K and 155 K at 300 K |38[]). The crystals in the DAMA experiment are at 
20 °C, i.e. T = 293.15 K; = 22.9 amu and Mj = 126.9 amu, thus M = 74.9 amu. The 
vibration amphtude ui we get using this value of M is plotted in Fig. ^ as a function of the 
temperature T. At room temperature (20 °C) it is ui = 0.0146 nm which is similar to the 

= 0.0145 nm 



while measured separate values of 



u^) for Na 



measured value of \/{u^ 
and I (always at room temperature) are 0.018 nm and 0.015 nm respectively. (To use the 
data in Ref. p0| we must take into account that the Debye- Waller factor B is B = 87r^(nf)). 




I 0.015 



0.010 



200 400 600 
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Figure 3: Plot of ui{T) for NaT (Eq. ||2|with M = (MNa + M/)/2) 



At low temperatures (T ^ 0) the individual vibration amplitudes become progressively 
different |41]. We are not taking this difference into account in our approach. However, at 
these temperatures the channeling fractions become so small (as we show below) that a better 
calculation is not important. 

In principle there are modifications to the continuum potentials due to thermal effects, 
but we are going to take into account thermal effects in the crystal through a modification 
of the critical distances which was found originally by Morgan and Van Vliet p7| and later 



by Hobler |15] to provide good agreement with simulations and data. For axial channels it 
consists of taking the temperature corrected critical distance rc(T) to be. 



r,(T) = y'rm + iciUiiTW, 



(3.29) 



where the dimensionless factor ci in different references is a number between 1 and 2 (see e.g. 



Eq. 2.32 of Ref. [|8| and Eq. 4.13 of the 1971 Ref .|231). 

For planar channels the situation is more complicated, because some references give a 
linear and other a quadratic relation between Xc{T) and ui. Following Hobler [^] we use an 



- 13 - 



<100> and |100| channels 
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Figure 4: (a) Static critical distances of approach and ui at 20 °C and (b) critical channeling angles 
at 20 °C with ci = C2 = 1 as a function of the energy of propagating Na (green/gray) and I (black) 
ions in the <100> axial and {100} planar channels. Here dach = c^pch- 



<111> and 1111} channels <111> and 1 1 1 1 ) channels 
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Figure 5: Same as Fig. ^ but for the <111> axial and {111} planar channels. Here dach 7^ dpch- 



equation similar to that for axial channels, 



(3.30) 



where again C2 is a number between 1 and 2 (for example Barret |14] finds C2 = 1.6 at high 
energies, and Hobler uses C2 = 2). We will mostly use ci = C2 = 1 in the following, to try 
to produce upper bounds on the channeling fractions. 

Using the T-corrected critical distances rc(T) and Xc{T) instead of the static lattice 
critical distances Tc and Xc in Eqs. and we obtain the T-corrected critical axial and 
planar angles. 

The static axial and planar critical distances are presented in Figs, ^(a) and ^(a) for the 
100 and 111 channels, respectively, together with the amplitude of thermal vibrations ui at 
20 °C. Figs. §(b) and |5|(b) show the temperature corrected axial and planar critical angles 
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<100> axial channel, Na recoil 



1 100) planar channel, Na recoil 
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E (keV) E (keV) 

Figure 6: Static and T corrected critical angles as a function of energy of Na recoil for T=77.2 K, 
T==20 °C, and T==600 °C for (a) <100> axial and (b) {100} planar channels. 



at 20 °C (with ci = C2 = 1) for the same channels as functions of energy of the traveling Na 
and I ions. We can clearly see in Fig. |5| that the critical angles become zero at low enough 
energies (for which the critical distance of approach should be larger than the radius of the 
channel) indicating the range of energies for which no channeling is possible. Figs, ^(a) and 
^(b) show the static and T-corrected critical angles at several temperatures for traveling Na 
ions in the 100 axial and planar channels respectively. 



4. Channeling of incident particles 

The channeling of ions in a crystal depends not only on the angle their initial trajectory 
makes with strings or planes in the crystal, but also on their initial position. Ions which start 
their motion close to the center of a channel, far from a string or plane, where they make an 
angle ip or ijjP respectively, defined in Eqs. p. 12 and 3.23, are channeled if the angle is smaller 



than a critical angle (as explain earlier) and are not channeled otherwise. Particles which 
start their motion in the middle of a channel (as opposed to a lattice site) must be incident 
upon the crystal (thus the title of this section). 

Here we show that to a good approximation we can use analytic calculations and repro- 



duce the channeling fraction in Nal presented in Ref. |12]. It must be noticed, however, that 
in Ref. the channeling fraction is computed as if the Na or I ions started their motion 
already within a channel (in fact close to the middle of the channel, where they assume the 
potential to be negligible), instead of starting from crystal lattice sites, as is the case in direct 
dark matter detection. Thus these calculations do not apply to direct dark matter detection 
experiments. 



Ref. |12] considers only a static lattice (i.e. no temperature effects taken into account) 



and the condition for axial channeling they use is ^ < ^2 [^i where is defined in Eq. |3.18| . 
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The equivalent condition for planar channeling used in Ref. |jT^ is 



il^P < 0pi = a^/Nd'p (ZiZaeV^a) 



1/3 



(4.1) 



The 9 pi angle is very similar to the planar critical angle found by Matyukhin [35| that we call 
■0?*^ and is given in Eq. C.4, Opi = (67r)~^/^'0c^''^. 



s 
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Figure 7: Maximum distance XniaxiE) traveled by channeled Na (green/gray) or I (black) ions in 
mixed channels of a Nal crystal (<100> and <111> axial and {100} and {110} planar channels). 

For an incident angle V with respect to each of the channels and an ion energy the 
fraction XinciE^il^) of channeled incident ions for axial and planar channels is Xinc = 1 if V' 
is smaller than the critical angle for the corresponding channel and zero otherwise. The Nal 
structure and all the different channels are explained in Appendix A. 

To find the total fraction Pjnc of channeled incident nuclei, we average Xinc over the 
incident direction q, 



Pu 



XmciE,q)dnq. 



(4.2) 



This integral cannot be solved analytically, so we integrated numerically by performing a 
Riemann sum once the sphere of directions has been divided using a Hierarchical Equal Area 
iso-Latitude Pixelization (HEALPix) fi^ (see Appendix B). 

A channeled ion can be pushed out of a channel by an interaction with an impurity such 
as the atoms of Tl in Nal (Tl). The probability density for an ion to find an impurity after 
propagating a distance x within the crystal is 



p[x) = iexp (^-j 



(4.3) 
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where A is the average distance between the Tl atoms. We take for A the value used by the 
DAMA collaboration, i.e. A = 120 nm |12|, which according to Ref. [|l^] corresponds to a 



molar concentration of 0.0013 Tl atoms for every Na atom. 

Here we will simply assume that if a channeled ion interacts with a Tl atom it becomes 
dechanneled and thus it does not contribute to the fully channeled fraction any longer. We 
thus neglect the possibility that after the interaction the ion may reenter into a channel, either 
the same or another, as we also neglect the possibility that initially non-channeled ions may 
be scattered into a channel. Both effects would increase somewhat the amount of channeled 
ions, but we do not have an analytic method of including them in our calculation. Thus the 
channeled fraction is simply reduced by the probability that the ion does not interact with a 
Tl atom. 



P'\E) = exp j^-^^i^J (4.4) 

Here Xmax(^) is the range of the propagating ion, i.e. the maximum distance a channeled ion 
with initial energy E can propagate along the channel. Within the channel the ion looses 
energy into electrons. We use the Lindhard-Scharff ^] model of electronic energy loss, 
valid for energies E < {Mi/2)zf^^VQ, where vq = e^/h = 2.2 x 10^ cm/sec is the Bohr's 
velocity |^]. Mi and Zi are the mass and charge of the propagating ion. This model is valid 
for E < 14.3 MeV for Na and E < 646.4 MeV for 1 in Nal. In this model the energy E{x) as a 
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Figure 8: (a) Channeling fraction for a 50 keV Na ion in different directions plotted on a sphere 
using the HEALPix pixelization; probability equal to one in red, and probability equal to zero in blue, 
(b) Fraction of channeled incident I (black) and Na (green/gray) ions as a function of their incident 
energy E with the static lattice without (dot dashed lines) and with (solid lines) dechanneling due to 
interactions with Tl impurities. The results of DAMA are also included (dashed lines). 

function of the propagated distance x and the initial energy E is the solution of the following 
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energy loss equation 



where v = y^2E/Mi is the ion velocity and K is the function 



= 3 — . (4.6) 



2 2 \ 2 

'3 I 73 



Here is a dimensionless constant of the order of Z^ [^] , A'^ is the number of atomic centers 
per unit volume and oq — 0.53 A is the Bohr radius of the hydrogen atom. Explicitly, an ion 
with initial energy E at x = has energy 

E{x) = e(i- (4.7) 



after traveling a distance x. The range of the propagating ion is 

Xr^UE) = (4.8) 

Fig. 1^ shows the maximum distance a^max traveled by channeled Na or I ions in mixed channels 
of an Nal crystal (<100> and <111> axial and {100} and {110} planar channels). The 
average distance A between Tl atoms is also shown. 

Fig. ^(a) shows the axial and planar channels of the Nal crystal in the HEALPix pix- 
elization of the sphere for incoming Na ions with an energy of 50 keV: red points indicate a 
channeling probability of 1 (when the incident angle is smaller than the critical angle with re- 
spect to any axial or planar channel) and blue points indicate a channeling probability of zero 
(when the incident angle is larger than the critical angle) . We include here only the channels 
with lower crystallographic indices, i.e. 100, 110 and 111, which provide the dominant con- 



tribution to the channeling fraction, as is also done in Ref. [12|. Fig. g(b) shows the fraction 
of incident I (black lines) or Na (green/gray lines) ions as a function of their incident energy 
E using the static lattice (solid lines). For comparison. Fig. ^(b) also shows the channeling 
fraction obtained by DAMA (dashed lines). Good agreement with the channeling fractions 
of DAMA is achieved only when dechanneling due to the interaction with Tl impurities is 
included. 



5. Channeling of recoiling lattice nuclei 

The recoiling nuclei start initially from lattice sites (or very close to them), thus blocking 
effects are important. In fact, as argued originally by Lindhard |^], in a perfect lattice and in 
the absence of energy-loss processes the probability that a particle starting from a lattice site 
is channeled would be zero. The argument uses statistical mechanics in which the probability 
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of particle paths related by time-reversal is the same. For example, in optics if a source of 
radiation and a point of observation are interchanged, the intensity of the light measured at 
the new place of observation is the same as the old. Thus the probability of an incoming 
ion to have a particular path within the crystal is the same as the probability of the same 
ion to move backwards along the same path 0. This is what Lindhard called the "Rule of 
Reversibility." 

Using this rule, since the probability of an incoming channeled ion to get very close to 
a lattice site is zero (otherwise it would suffer a large angle scattering and it would not be 
channeled), the probability of the same ion to move in the time-reversed path, starting at a 
nuclear site and ending inside a channel, is zero too. However, any departure of the actual 
lattice from a perfect lattice, for example due to vibrations of the atoms in the lattice, would 
violate the conditions of this argument and allow for some of the recoiling lattice nuclei to be 
channeled. 

The channeling of particles emitted at lattice sites due to lattice vibrations, such as 
protons scattered at large angles, was measured and already understood in the 70's. Komaki 
et al. g in a 1971 paper titled "Channeling Effects in the Blocking Phenomena" observed 



channeling of protons scattered at large angles within thin Si and Ge crystals and explained 
it as due to the fact that the scattering or emitting lattice atom is not exactly at the lattice 
site because of thermal vibrations." They fit their data using the model presented by Komaki 



and Fujimoto |3C] one year earlier. 



We now estimate the channeling fractions using the formalism presented so far. 
5.1 Channeling fraction for each channel 

We need to know how probable it is for the recoiling nucleus to be at a particular distance r 
from its equilibrium position in a crystal row when it collides with a WIMP. The probability 
distribution function g{r) of the perpendicular distance to the row of the colliding atom due to 
thermal vibrations can be represented by a two-dimensional Gaussian (as done by Lindhard 
and many others Q, the relevant vibrations being in the plane orthogonal to the row), 

<7(r) = 4exp(-rV2u?). (5.1) 



The one dimensional vibration amplitude ui is given in Eq. |3.26 . 

The channeled fraction Xa.na.i{E, q) of nuclei with recoil energy E moving initially in the 
direction q making an angle (j) with respect to the axis is given by the fraction of nuclei 
which can be found at a distance r larger than a minimum distance rj^min from the row at 
the moment of collision, determined by the critical distance of approach as shown in the next 
subsection, 

/oo 
drg{r) = exp {-rl^iJ2ul). (5.2) 
- i,min 

Note that here we are approximating the upper limit of the integral of g{r) with oo, instead 
of the radius of the axial channel (iach/2. This is a good approximation because (iach/2 — 10a 
or more and the integral is dominated by the values of g{r) close to ui <^ 2a. 
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If (/) > -0 c no channeling can occur and Xa.xiaiiE , (p) — 0. This can easily be seen from 
Eqs. 3.8, 3.11, 3.12 and \i.l6 , taking into account that U{ri + dtancj)) > U{rc). 

For a planar channel, the Gaussian thermal distribution for the planar potential is one- 
dimensional (the relevant vibrations occurring perpendicularly to the plane), 

g{x) = (27ruf)-^/2gxp(-xV2n?). (5.3) 

This is normalized to 1 for —oo < x < +oo. In our calculations we only consider positive 
values of x for each plane, thus we multiply g{x) by a factor of 2 to find the fraction of 
channeled nuclei for a planar channel, 

Xpianar(£^,<^) = / '^9{x)dx = —= y= dx = evfc [ ) . (5.4) 

Here (p is the angle q makes with the plane, defined as the complementary angle to the angle 
between q and the normal to the plane, or as the smallest angle between q and vectors lying 
on the plane. Similar to the axial case, we approximate the upper limit of the integral of g{x) 
with oo, instead of the radius of the planar channel (ipch/2. This is a good approximation 
because (ipch/2 ^ lOui or more, and erfc[(ipch/(2\/2ni)] is negligible. Also in this case, 
XpianariE, (/)) = if (/) is larger than the critical channeling angle of the particular channel, i.e. 
if (/> > V'?- 

We conclude this subsection by noticing an important point. In Eq. 5.2, rj min, which 
is a function of rc{T), enters exponentially. Thus any uncertainty in our modeling of rc{T) 
becomes exponentially enhanced in the channeling fraction. The same happens with the 
dependence of the channeling fraction in Eq. on Xi^mm, which depends on Xc{T). This is 
the major difficulty of the analytical approach we are following. 

5.2 Minimum initial distance of the recoiling lattice nucleus 



For axial channels, using Eqs. |3.8| , |3.11| , |3.12| and |3.1q , we can write the condition for chan- 
neling as 

Ssin2 + C/(ri + dtan0) = [/(r^in) < U{rc{E)). (5.5) 
Therefore, the minimum initial distance 

^i,min is the solutioii of the ec[ua.tioii 

f^(ri,min + c?tan</)) = U{rc{E)) - Esm^ (p. (5.6) 
Inverting the function U{r) we obtain 

r^,min(£^,</') +c^tan</. = U-\U{rc{E)) - Esm^cj)]. (5.7) 
The inverse of Lindhard's potential function U (r) is 

U-\r)= , (5.8) 
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which together with the expression for rc{E) in Eq. 3.14 yields a fuhy analytic expression for 



^iE,< 



Ca 



1 + ^) exp(-2sin2 0/^2)_i 



dtan (h. 



(5.9) 
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Figure 9: Upper bounds to the channeling fractions of Na recoils for single planar (green/gray lines) 
and axial (black lines) channels, as function of the recoil energy E, for T= 293 K and ci = C2 = 1, 
without including dechanneling. Two additional channels not included in the total channeling fractions 
are also shown: axial [211] (brown line) and planar (210) (cyan line). 




Applying the same arguments to planar channels, we have 

U{xi^jjiin + dp tan cf)) = U {xc{E)) — E svo? (j) 
For Lindhard's potential, the minimum initial distance is given by 



(5.10) 



.iE,< 



a 
2 



+ C2-^-sin2</,/V;2 



ft + C72-f -sin^cA/V^ 



dp tan ( 



(5.11) 



Here, Xc{E) is found in Eq. 3.22. 

Fig. |9| shows upper bounds to the channehng fractions of Na recoils for individual channels 
with ci = C2 = 1 and T= 293 K, without including dechannehng. The black and green 
(or gray) lines correspond to single axial and planar channels respectively. Two additional 
channels not included in the total channeling fractions are also shown for comparison: the axial 
[211] (brown line) and the planar (210) (cyan line) channels. The upper bounds of channeling 
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Figure 10: Channeling probability Xrcc{E,ci) (Eq. 5.13 ) for a 200 keV recoil of (a) an Na ion and 
(b) an I ion at 20 °C with ci = C2 = 1 and neglecting dechannneling. The probability is computed for 
each direction and plotted on a sphere using the HEALPix pixelization. The red, pink, dark blue and 
light blue colors indicate a channeling probability of 1, 0.625, 0.25 and zero, respectively. 



fractions for planar channels are more generous than those of axial channels because of our 



choice of Xc in Eq. 3.22. This does not mean that planar channels are dominant in the actual 
channeling fractions. 

5.3 Channeling fraction 

The geometric channeling fraction is the fraction of recoiling ions that propagate in the 1st, 
or 2nd, or ... or 26th channel. Here "geometric" refers to assuming that the distribution 
of recoil directions is isotropic. In reality, in a dark matter direct detection experiment, 
the distribution of recoil directions is expected to be peaked in the direction of the average 
WIMP flow. For comparison with previous work of others, here we examine this geometrical 
channeling fraction, and postpone the case of a WIMP wind to another paper [p5|]. 

We include only the channels with lowest crystallographic indices, i.e. 100, 110 and 111, 
which are in total 26 axial and planar channels, as explained in Appendix A. We have also 
checked other axial and planar channels, such as the [211] and (210) channels shown in Fig. ^, 
and found that their contribution to the channeling fractions is negligible (additional planar 
channels are always less important than the planar channels we keep and the same happens 
for axial channels). 

The probability XTeciE,ci) that an ion with initial energy E is channeled in a given 
direction q is the probability that the recoiling ion enters any of the available channels. We 
compute it using a recursion of the addition rule in probability theory over all axial and 
planar channels: 

P(Ai or A2) = P(Ai) + P(A2) - P{Ai)P{A2). 
P{Ai or A2 or A3) = P{Ai or A2) + ^(^3) - P{Ai or ^2)^(^3)- (5.12) 
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We continue this recursive computation until we find the probabihty with which the recoihng 
ion goes into any of the 26 channels 



Xrec{E, q) = P{Ai 0VA2 or ... or A 



26; 



(5.13) 



For each channel {k = 1, . . . ,26), P{Ak) = Xaxiai-fc(-E^, 0fc) or P{Ak) = Xp\a.ni,r-k{EAk) 
for an axial or planar channel, respectively. Notice that P{Ak) 7^ only for the channels for 
which < {ipc)ki for which the angle that q makes with the axis or plane of the channel, 
respectively, is smaller than the critical angle for the channel. 

Here we are treating channeling along different channels as independent events, so that 
the conditional probabilities coincide with the non-conditional probabilities, e.g. P(j4i|A2) = 
P{Ai). This is correct for different axial channels, which never overlap, and a good approxi- 
mation for different planar channels. However, axial channels happen at the crossing of two 
or more planar channels, thus channeling into axial and planar channels may not be entirely 
independent. We prove in Appendix D that considering them as independent is, however, a 
good approximation. 
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Figure 11: Upper bounds to the channeling fraction of Na and I recoils as a function of the recoil 
energy E for T=600 °C (green/light gray), 293 K (black), and 77.2 K (orange/dark gray) in the 



approximation of ci = C2 = 1, (a) without and (b) with dechanneling as in Eq. 4.4 



Fig. 10 shows the channeling probability for an E = 200 keV recoil of (a) Na or (b) I 
at 20 °C with ci = C2 = 1 and neglecting dechanneling, computed for each direction q and 
plotted on a sphere using the HEALPix pixelization. The red, pink, dark blue and light blue 
colors indicate a channeling probability of 1, 0.625, 0.25 and zero, respectively. 

To obtain the geometrical channeling fraction, we average the channeling probability 
XTec{E,c[) over the directions q, assuming an isotropic distribution of the initial recoiling 
directions q, 

PrUE) = ^ / Xrec{E,Cl)dng. (5.14) 



This integral is computed using HEALPix [^] (see Appendix B). 
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Figure 12: Same as Fig. |ll| but for ci — C2 ~ 2. 
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Figure 13: Same as Fig. ^ but for ci = C2 = (static lattice), provided as an upper bound with 
respect to any non-zero values of ci and C2. 



The channeling fraction as a function of recoil energy is shown in Figs. 11 and 12 with 
ci = C2 = 1 and ci = C2 = 2, respectively. These curves include thermal effects in the 
lattice at various crystal temperatures, and are shown (a) without and (b) with dechanneling 
according to Eq. Note that the dechanneling as included here is possibly too extreme, 
since it does not allow for the possibility of an ion reentering a channel (the original or 
another one) after a collision with a Tl impurity. As shown in Fig. pT|(a), the maximum of 
the channeling probability for Na recoils is at energies close to 200 keV. 

Notice that while the effect of increasing temperatures on the initial position of the 
recoiling nucleus makes the channeling fractions larger (the recoiling nucleus can be initially 
further out from the string, i.e. ui in Eqs. 5.2 and 5.4 increases), the effect of increasing 
temperatures on the other lattice atoms is to increase the critical distances, which makes 
the channeling fractions smaller (rj^min and Xj^min in Eqs. |5.2| and 5A increase). Figs. 11 and 
|l^ show that in our calculations the first effect is almost always dominant, except that for 
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ci = C2 = 2 at some energies the temperature effects in the lattice are larger (see the left 



panel of Fig. 12, where some fractions are smaller at higher temperatures). Neglecting the 
temperature effects in the lattice, by setting ci = C2 = 0, but including the thermal vibrations 
of the nucleus that is going to recoil, we obtain the largest estimates for the channeling 
fractions. Although it is physically inconsistent to take only the temperature effects on the 
initial position of the recoiling nuclei but not on the lattice, this was done by Lindhard Q and 
Andersen early on, and we do it here because it provides the most generous upper bound 
on the channeling fraction (any non-zero value of ci or C2 would lead to smaller fractions). 



6. Conclusions 



We have studied the channeling of ions recoiling after collisions with WIMPs inside an Nal 
(Tl) crystal. Channeled ions move within the crystal along symmetry axes and planes and 
suffer a series of small-angle scatterings that maintain them in the open "channels" between 
the rows or planes of lattice atoms and thus penetrate much further into the crystal than in 
other directions. In order for the scattering to happen at small enough angles, the propagating 
ion must not approach a row or plane closer than a critical distance Tc for axial or Xc for planar 
channels. For a "static lattice" that here means a perfect lattice in which all vibrations are 



neglected, rc and Xc are given in Eqs. 3.14 and 3.22 



The channeling of ions in a crystal depends not only on the angle their initial trajectory 
makes with rows or planes in the crystal, but also on their initial position. Ions which start 
their motion close to the center of a channel, far from a row or plane, at an initial angle ip or 
ipP (see Eqs. 3.12| and p. 23 ), are channeled if the initial angle is smaller than the critical angle 



in Eqs. 3.17 and 3.25, respectively, and are not channeled otherwise. We have found that the 



channeling of lattice ions recoiling after a collision with a WIMP is very different from the 
channeling of incident ions, and that the fraction of recoiling lattice ions that are channeled 
is smaller. 

The nuclei ejected from their lattice sites by WIMP collisions are initially part of a row or 
plane. They start from lattice sites or very close to them, thus blocking effects are important. 
In fact, as argued originally by Lindhard [^, in a perfect lattice and in the absence of energy- 
loss processes, the probability that a particle starting from a lattice site is channeled would 
be zero. This is what Lindhard called the "Rule of Reversibility." However, any departure 
of the actual lattice from a perfect lattice due to vibrations of the atom, which are always 
present, violate the conditions of this argument and allow for some of the recoiling lattice 
nuclei to be channeled. Thus, the channeling fraction of recoiling ions is very temperature 
dependent. 

Due to vibrations in the crystal, the atom that interacts with a WIMP may be displaced 



from its position in a perfect lattice with a probability given in Eqs. 5.1 and 5.2. It is 



this displacement which allows for a non-zero probability of channeling, given in Eqs. 5.2 



and 5.4. At high temperatures, the atoms vibrate with larger amplitudes, the recoiling ion 



can start further away from a row or plane (i.e ui in Eqs. 5^ and ^.4| is larger), and the 
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channeling fractions increase. This is illustrated in Figs. 11, 12 and 13. These three figures 
differ because the lattice vibrations (of all the other atoms in the crystal, besides the recoiling 
one) increase the critical distances of approach and reduce the critical angles for channeling 
(unless ci = C2 = 0), which in turn decreases the channeling fractions as the temperature 
increases (rj^min and Xi,min in Eqs. 5.2 and ^jl increase). 



Without better data or simulations of Na and I ions propagating in a Nal crystal, the 
best we can do is to write the T-corrected minimal distances of approach rc{T) and Xc{T) as 
functions of the two parameters ci and C2- The values of ci and C2 we found in the literature 
for different materials and propagating ions are between 1 and 2 (see Eq. ^and^). 
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Figure 14: Channeling fractions at T=293 K for Na (solid lines) and 1 (dashed lines) ions for 
c = ci = C2 — 1 (black) and c ~ ci = C2 = 2 (green/gray) cases (a) without and (b) with dechanneling 



included as in Eq. 4.4 



The coefficients ci and C2 determine the dependence of the critical distances of approach 
on the vibration amplitude ui{T). This dependence is important at relatively large energies, 
when the static critical distances Tc and Xc become small with respect to ui. 

We are already providing upper limits to the channeling fraction due to our choice of Xc 



(see discussion after Eq. 3.22 ). In Fig. 13 we decided to show absolute upper bounds to this 
probability in the unrealistic case in which temperature effects are only taken into account 
in the vibrations of the atom interacting with the WIMP, but not on the other atoms in the 
lattice (so ci = C2 = 0). More realistic upper bounds to the channeling fraction (without 
dechanneling) are given in Fig. 11, in which ci = C2 = 1. The 20 °C curve from this figure 
is displayed again in Fig. 14, in which we show what we consider to be our best results. The 
larger are ci and C2, the larger are the temperature corrected critical channeling distances, the 
smaller are the temperature corrected critical channeling angles, and thus the smaller are the 
resulting channeling fractions. If ci = C2 = 2, found by Hobler to be necessary to reproduce 
the measured critical angles at < 100 keV in Si (see our companion paper devoted to Si 



and Ge [46|), the channeling fractions are smaller, as shown in Fig. 12 from which the 20 °C 



curve is copied in Fig. 14 
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Fig. 14(a) shows what we consider to be our main predictions for the range expected as an 
upper Hmit to the channehng fraction in Nal (Tl), if dechannehng is ignored. DechanneHng 
happens when the channeled ion encounters impurities or defects. Fig. 14(b) shows the 
channehng fraction reduced by the probabihty of the channehng ion to not interact with a Tl 
atom (see Eq. 4.4). This way of taking into account dechannehng may be too extreme, as it 
neglects the probability that the ion after the collision with a Tl atom may reenter a channel 
(either the same channel or another one) and be again channeled. 



As we see in Fig. 14(a), neglecting interactions with Tl atoms, the channeling fraction 
is never larger than 5% and the maximum happens at lOO's of keV. This maximum occurs 
because the critical distances decrease with the ion energy E, making channeling more prob- 
able, and the critical angles also decrease with E, making channeling less probable. All this 
is without dechanneling. 

With dechanneling, the probability that the channeled ion does not interact with a Tl 
atom decreases with energy (since more energetic ions propagate further within channels). 
Thus, interactions with Tl atoms decrease the channeling fraction at high energies. The 
simple extreme model of dechanneling we use in this paper predicts much smaller fractions, 
at most in the 0.1% level, with the maximum shifted to small energies, less than 10 keV (see 
Fig. ^(b)). This reduction may eventually prove to be too extreme and at present we do not 
have a better formalism to model dechanneling. 

With the simple model of dechanneling we used we could reproduce the channeling frac- 
tions computed by the DAMA collaboration (which, however, apply to ions which start their 
motion close to the middle of a channel and not to the case of WIMP direct detection). 
For completeness, in Appendix C we also include a model of planar channeling proposed by 
Matyukhin |35| , which produces larger channeling probabilities, but which we think is not 
correct. 

The analytical approach used here can successfully describe qualitative features of the 
channeling and blocking effects, but should be complemented by data fitting of parameters 
and by simulations to obtain a good quantitative description too. Thus our results should in 
the last instance be checked by using some of the many sophisticated Monte Carlo simulation 
programs implementing the binary collision approach or mixed approaches. 
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A. Crystal structure and other data for Nal 



Nal is a diatomic compound that has two interpenetrating face-centered cubic (f.c.c.) lattice 
structures displaced by half of a lattice constant with 8 atoms per unit cell. The lattice 
constant aiat of a cubic crystal system refers to the constant distance between unit cells of 
one of the f.c.c. lattices in the crystal, and for Nal it is aiat =0.6473 nm at room temperature 



(Table 3.4 in Appleton and Foti [^2[)- -^iS- shows one eights of the unit cell of the Nal 
crystal. The red and blue spheres represent Na and I ions respectively. The shortest distance 
between Na and I ions in Fig. |l^ is half the lattice constant, aiat/2. 




Figure 15: One eights of the Nal crystal unit cell with the red and blue spheres representing Na and 
I ions respectively. The sohd, dashed and dot-dashed lines show the <100>, <110>, and <111> axes 
respectively. The {100}, {110} and {111} planes are perpendicular to the respective axes with equal 
indices. 



The atomic mass and atomic number of Na and I are = 22.9 amu, Mj = 126.9 amu, 
Znsl = 11 and Z[ = 53. When computing -01 in Eq. ^]2| for Na recoils, we take Zi = Z^a and 
Z2 equal to an effective atomic number of the row or plane in the channel, which depends on 
the composition of the row or plane. "Mixed" channels, for example the rows < 100 > and 
< 111 >, or the planes {100} and {110}, contain both Na and I ions in alternation; they have 
Z2 = Z = (.^Na + -^i)/2. "Pure" channels, for example the row < 110 > or the plane {111}, 
contain atoms of a single species, only Na or only I; they have Z2 = .^Na or Z2 = Zi. Thus, 
for Na recoils from the row where it originally was, we have 

sin^<^ = ^^^|^. (A.1) 
Similarly, for I recoils Zi = Zi, and for mixed channels Z2 = Z = {Z^i^ + Zi) /2 while for pure 
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channels we use Z2 = Z\. We have 



-^W = ^^. (A.2) 



With respect to the Thomas-Fermi screening distance, for Na recoils from a mixed row 
or plane we use the average 

ONa = (oNaNa + aNai)/2 = 0.01149 nm, (A. 3) 

where ONaNa = 0.4685(zi{^ + z]IIY'^I^ = 0.01327 nm and ONai = 0.4685(^^1^ + zl'\^l^ = 
0.009711 nm correspond to an Na scattering from an Na and an I lattice atom, respectively. 
On the other hand, for Na recoils from a pure row or plane we use aNaNa because the row or 
plane from which the recoiling Na ion was emitted contains only Na atoms. Similarly, for I 
recoils from a mixed row or plane, we use 

ai = (an + aNai)/2 = 0.008784 nm, (A.4) 

where an = 0.4685(z//V^J^^)"^/^ = 0.007857 nm and CNal correspond to an I ion scattering 
from an I and an Na lattice atom, respectively. For I recoils from a pure row or plane we use 
an , since the row or plane the recoiling ion is emitted from is made of I ions only. 
Finally, to compute ipa for Na recoils 

sm^.= (^^^!^%^)\ (A.5) 

where Z2 = Z and a = ONa or Z2 = ^Na and a = ONaNa for mixed or pure rows and planes 
respectively. For I recoils 

smV'a = , (A.6) 



where Z2 = Z and a = ai or Z2 = Z^^ and a = an for mixed or pure rows and planes 
respectively. 

We here review the crystallographic notation for directions in the lattice. Once an origin 
of the coordinate system is fixed on a lattice point O, any position vector of a point on the 
crystal lattice can be written as R = nia + n2h + n-^c with ni, ^2, and ns specific integer 
numbers. The vectors a, b, and c are the basis vectors of the crystal lattice, and are three 
noncoplanar vectors joining the lattice point O to its near neighbors. For the cubic lattice of 
Nal, the three vectors a, b, c form a Cartesian frame and their length is aiat/2 (they are the 
sides of the cube in Fig. 15). The integers ni, n2, and 713 can be positive, negative, or zero. 



The direction of a crystal axis pointing in the direction R is specified by the triplet [nin2?T-3] 
written in square brackets, when ni, 712, and 71.3 are positive or zero. Note that if there is a 
common factor in the numbers ni, 722, 71.3, this factor is removed. Moreover, negative integers 
are denoted with a bar over the number, e.g. —1 is denoted as I and the —y axis is [010] 
direction. Fig. |l^ shows the directions of the [100], [110] and [111] axes. 
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In a cubic crystal, because of the symmetry of the unit cell, the directions [100], [010], 
and [001] are equivalent. All directions equivalent to the [nin2n3] direction are denoted by 
<nin2n3> in angular brackets. For example, <100> indicates all six directions [100], [010], 
[001], [100], [010], and [001]. The plane perpendicular to the [^1722/13] axis is denoted by 
(nin2n3). For example, the plane perpendicular to the [100] axis is denoted by (100), and 
that perpendicular to [101] by (101). The integers ni, 712, and are called Miller indices. 

When the unit cell has cubic symmetry, we can indicate all planes that are equivalent 
to the plane (hkl) by curly brackets {hkl}. For example, the indices {100} refer to the six 
planes (100), (010), (001), (100), (010), and (001). The negative sign over a number denotes 
that the plane cuts the axis on the negative side of the origin. 

We will only consider the lower index crystallographic axis and planes. For axial chan- 
neling we will consider the <100>, <110> and <111> axes and for planar channeling we 
consider the {100}, {110} and {111} planes perpendicular to them. 

To compute the interatomic spacing d in axial directions and the interplanar spacing dp^h 
in planar directions, we have to multiply the lattice constant by the following coefficients [0|: 

• Axis: < 100 >: 1/2 , < 110 >: 1/^2 , < 111 >: y/3/2 

• Plane: {100} : 1/2 , {110} : 1/2^2 , {111} : 1/2^3 

For Nal, the Debye temperature is = 165 K, and the crystals in the DAMA experiment 
are at a temperature of 20 °C or 293.15 K. 

B. HEALPix Pixelization 

The Hierarchical Equal Area iso-Latitude Pixelization (HEALPix) [^2| provides a convenient 
way of dividing the surface of a sphere into equal area sectors. An integral over directions can 
then be performed as a simple Riemann sum. HEALPix has been introduced to pixelize data 
on a sphere and has been used by cosmic microwave background experiments like WMAP 
and BOOMERANG. 

In HEALPix, the base resolution comprises 12 pixels in three rings: one ring around the 
north cap, one ring around the south cap, and one ring around the equator. At a higher 
resolution, each base pixel in each ring is divided into smaller pixels of equal area. The 
resolution parameter of the grid is Ngidc and it defines the number of divisions along the side 
of a base-resolution pixel which is needed to obtain a partition with higher resolution. We 
choose the resolution parameter of the grid to be A^side = 50. 

A HEALPix map has A'pixei = 12N^^^^ pixels, each with the same area. The angular 
resolution of the map can be estimated by computing the solid angle covered by each pixel 
= 47r/A'^pixei, and finding the typical diameter of each pixel as if it were small and of circular 
shape, 6*1.08 = 2y/Q/-K. This gives 0res = 2/(V3iVside) = 66.2°/iVsido- By choosing iVsido = 50, 
we have 30,000 pixels on the sphere, and a resolution of 1.3 degrees. If the HEALPix is 
properly aligned with the cubic crystal so that there is a pixel in each <100> direction, this 
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resolution should be sufficient for computing our integrals accurately. We have tried different 
values of N^ide, up to N^ide = 400 (for which the resolution is ^res — 0.166°), and we found 
that the value of the channeling fraction already converges within one percent for Ngide = 20. 
Thus, we used A^side = 50 as a safe value in out calculations. 

The following algorithm was used to generate a list of unit vectors on a sphere, with 
each unit vector in the direction of one of the HEALPix pixels. Let p be the pixel index, with 
p = 0,1, . . . , A'pixei — 1- We start with the definitions: 
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We proceed to compute the cylindrical coordinates z and (p of the direction of the p-th 
HEALPix pixel. If g > p*^^^""^^^^, the pixel belongs to one of the polar rings and we succes- 
sively compute 
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Here [xj is the minimum integer less or equal to x, and x mod y is the remainder of the 
integer division of x by y. li q < the pixel belongs to the equatorial ring and we 
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Finally, the direction vector of the p-th pixel is given by 



lip = — z'^ cos if, \J\ — sine/?, z\ . (B.16) 



C. Matyukhin model of planar channeling 

Matyukhin in Ref [^] uses a different condition for planar channeling than we used previously 



which consists of Eq. 3^ written in terms of the planar potential, namely Up{x) < 8E/dp 
(where " denotes the second derivative with respect to x). We have not been able either 
to derive this condition for the planar potential or to find the derivation of this condition 
anywhere. We have not found the results derived from this condition compared with data 
either, but we mention the model for completeness. For Lindhard's planar potential, the 
condition used by Matyukhin becomes 

8 [^l^^ + C^a^]V2^ ^^■') 



where Xmin is the minimum distance of approach to the plane, dp = and 



E2 = Eipl = 2T:NdpZiZ2e^a. (C.2) 



From Eq. |C.1| , the smallest possible value of Xmin is now 

For all the energies we consider here (~ keV and above), x^^ is smaller than dpch/2, thus 
Up{dpch/'i) can be neglected in the definition of the critical planar angle V'c , Eq. 3.25| . At low 



energies^ < dp^ E2{SC a?)-^ wehavexf (^) ~ Ca{dp^E2/SCa'^EYl'^ = {C'^T:ZiZ2e^a? /AEY'^ 
and 



. dp E J \ E 

where tl;f \E) = (67r)i/3^p/ and 9pi is the critical angle used by DAMA. 

As E approaches the value dp'^ E2{^C a?)~^ , x^^ {E) approaches zero. At larger energies, 
E > dp^ E2{^Ca?)~^ , x^^ {E) in Eq. p.3| becomes imaginary, but x^^ {E) is by definition real 
and positive thus one could take it to be zero. Matyukhin takes in this case x^'\E) = a 
instead. Either way, in this energy range ipc\E) ~ il^a which is the value given by Lindhard 
for the "breakthrough" angle [0] necessary to have -Eperp = Up{0). We take x^^ (E) = a 
wherever the prediction of Eq. p.3| is smaller than a. Including temperature corrections due 
to the thermal and zero point energy vibrations of the atoms in the lattice, we have Xc{T) as 
in Eq. g. 



The equations of Matyukhin coincide with those presented here if C = 1 (but, following 
Lindhard, we take C = \/3 instead). 
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Figure 16: Comparison of |T^(a) (top left) critical distances of approacli and ^(b) (top riglit) critical 
angles at 20 °C witli ci = C2 = 1 in tlie {100} planar channel predicted by our main model (solid lines, 
see Figs. ^, ||) and by Matyukhin's (dashed lines). |l^(c) (bottom left) and |l^(d) (bottom right), same 
for the {111} planar channel. Green/gray lines are for Na and black for I propagating ions. 



A comparison of the static critical distances of approach Xc{E) in our method (using 
Eq. 3.22 and 3.25 ) and in Matyukhin's model, and of the critical angles for ci = C2 = 1 
in both models is shown in Figs. 16(a) and 16(c) and in Figs. |l^(b) and |l6|(d) respectively 
for the {100} and {111} planar channels respectively. The Matyukhin critical distances of 
approach are smaller (and thus the critical angles larger) than those in our main method 
at low energies, which leads to higher channeling fractions, as shown in Figs. ^ and 18 for 
ci = C2 = 1 and for ci = C2 = 2 respectively. In these figures the left panels are without 
and the right panels with dechanneling as in Eq. 4.4 included. The channeling fractions using 
Matyukhin's model are much higher than the fractions we obtain with our method, but, as 
explained above, we do not trust Matyukhin's model. The critical distances in Matyukhin's 
model (see Figs, ^(a) and 0(c)) have a discontinuous slope at the energy where they become 
constant and this shows in the channeling fraction curves also as discontinuities in slope (at 
the values of E at which different important channels have a sharp change in x^^). 
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Figure 17: Channeling fractions using Matyukhin's model for the planar channel as a function of 
the energy of recoiling Na and 1 ions for T=600 °C (green/light gray), 293 K (black) and 77.2 K 
(orange/dark gray) for T-corrections included in the lattice with c\ — ci ~ \^ (a) without and (b) 



with dechanneling included as in Eq. 1.4. The probabilities are larger than in our main method, but 
we do not trust Matyukhin's approach. 
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Figure 18: Same as Fig. |l^ (for Matyukhin's model) but for ci = C2 = 2. 
D. Probability of correlated channels 



In Eq. |l2 we treat channeling along different channels as independent events. Here we prove 
that this procedure is adequate for our purpose of providing upper bounds to the channeling 
fractions. 



For an axial channel when (j) < ipc (otherwise Xaxiai = 0), the integration region in Eq. 5.2 
is the exterior of an infinitely long cylinder of radius rj^niin(-E') 4') and axis coincident with 
the channel axis. Similarly, for a planar channel when (j) < if)^ (otherwise Xpianar = 0), the 
integration region in Eq. 5.4 is the exterior of an infinite slab of half-thickness Xi-camiE^cj)). 
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Let us consider the complements of the integration regions, i.e. the regions excluded in the 
integrals. These are the regions interior to a cylinder (for an axial channel) or a slab (for a 
planar channel). 

Note that only channels making angles with the direction q (of the initial momentum) 
smaller than their respective critical angles contribute to the union of integration regions. 
Therefore the problem of combining channels arises only when a recoil direction belongs to 
more than one channel, and this happens if the channels overlap for some directions of q. 
For the cases we consider (cubic lattices), only axial channels overlap with a subset of planar 
channels, or two or more planar channels overlap with each other. Notice that two different 
planar channels crossing at an angle, overlap in a parallelepiped of very long length on one 
side, thus one can define an inscribed cylinder within the parallelepiped, whose diameter is 
equal to the smallest of the two widths of both planar channels. 

We can obtain an upper limit to the channeling probability of overlapping channels by 
replacing the intersection of the complements of the integration regions with the inscribed 
circle of radius tmin equal to the minimum of the rj^min or Xj^min among the overlapping 

is 



channels. Then an upper bound to the probability Xrec(-E',q) in Eq. 5.14 

Xrec(-E,q) < / 



drg{r) = exp {-ry^i^/2u-^). 



(D.l) 



'^'MIN 



When only one channel is open (i.e. (j) < ipc for only one channel), we still use Eq. |5.2| or 5.4 
for the channeling probabilities. 

Fig. ^ shows the comparison of this upper limit with the fractions we computed using 



Eq. |5.12| . The two are practically indistinguishable. This proves that the method we used 
(see Section 5.3) is adequate for our purpose of providing upper bounds to the channeling 
fractions. 
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Figure 19: Maximum chamieling fractions using Eq. D.l for a — C2 — c and c — (dashed green), 
c = 1 (dashed yellow) and c — 2 (dashed cyan) compared with the results of our method of Section 
5.3 (solid black lines) for the same models for (a) Na ions and (b) 1 ions propagating in Nal crystal at 
T — 293 K. Notice that the lines overlap. 



- 35 - 



References 

[1] C. Savage, G. Gelmini, P. Gondolo and K. Freese, arXiv:1006.0972v2 [astro-ph.CO]. 

[2] F. T. Avignone, R. J. Creswick and S. Nussinov, arXiv:0807.3758 [hep-ph]; R. J. Creswick, 
S. Nussinov and F. T. Avignone, arXiv:1007.0214v2 [astro-ph.IM]. 

[3] P. K. Rol et al. "Proc. Fourth Int. Conf on Ionization Phenomena in Gases" ed. by N.R. 
Nilsson, North Holland, p.257 (1960). 

[4] M. T. Robinson Appl. Phys. Lett. 1, 49 (1962). 

[5] J. Stark, Phys. Z. 13, 913 (1912). 

[6] J. Lindhard, Kongel. Dan. Vidcnsk. Selsk., Mat.-Fys. Medd. 34, No. 14 (1965). 
[7] D. S. GemmcU, Rev. Mod. Phys. 46, 129 (1974). 

[8] H. Erramh and G. Blondiaux, Appl. Radiat. Isot. 46, 413 (1995); J. U Andersen et al., Phys. 
Rev. C 78, 064609 (2008). 

[9] C. Cohen and D. Dauvergne Nucl. Instrum. Methods Phys. Res. B 225, 40 (2004). 

[10] R. Altman, H. B. Dietrich, R. B. Murray, and T. J. Rock, Phys. Rev. B 7, 1743 (1973). 

[11] E. M. Drobyshcvski, Mod. Phys. Lett. A 23, 3077 (2008) [arXiv:0706.3095 [physics.ins-det]]. 

[12] R. Bernabei et al, Eur. Phys. J. C 53, 205 (2008) [arXiv:0710.0288 [astro-ph]]. 

[13] L. Rubin and J. Poate, The Industrial Physicist, June/July 2003, p. 12-15. 

[14] J. H. Barrett, Phys. Rev. B3, 1527 (1971). 

[15] MDRANGE, http://beani.acclab.helsinki.fi/knordlun/mdh/mdh_program.htnil; SRIM, http: 
//www.srim.org/; TRIM, J. F. Ziegler, "Ion Implantation Technology", Ion Implantation 
Technology Co. (1996); MARLOWE and UT-MARLOWE, Y. Chen et al, IEEE Trans. 
Electron Devices, vol. 49, no. 9, 1519 (2002); Crystal-TRIM, http://www.fzd.de/pls/rois/; 
REED, K. M. Beardmore and N Gronbech- Jensen, Phys. Rev. B 60, 12610 (1999). 

[16] V. Bykov et al. Nucl. Instrum. Methods Phys. Research (NIM)B 114, 371 (1996). 

[17] K. Nordlund, Comput. Mater. Sci. 3, 448 (1995). 

[18] G. Hobler, Radiation effects and defects in solids 139, 21 (1996); G. Hobler, Nucl. Instrum. 
Methods Phys. Research (NIM)B 115, 323 (1996). 

[19] K. Cho et al, Phys. Research B7/8 , 265 (1985). 

[20] C. J. Andreen and R.L. Hines, Phys. Rev. 159, 285 (1967). 

[21] I. Bergstrom et al.. Can. J. Phys. 46, 2679 (1968). 

[22] A. van Wijngaarten, B. Miremadi, N.M.A. Finney and J.N. Bradford, Phys. Rev.185, 490 

(1969). 

[23] K. M. Lui, Y. Kim, W. M. Lau and J. W. Rabalais, Jour. Apphed Phys 86, 5256 (1999). 
[24] S.M. Hogg, B. Pipclcers, A. Vantomme and M. Swart, Applied Phys. Lett. 80, 4363 (2002). 
[25] Z. L. Fang, W. M. Lau and J. W. Rabalais, Surface Science 581, 1 (2005). 



- 36 - 



[26] J. U. Andersen, Kongel. Dan. Vidensk. Selsk., Mat.-Fys. Medd. 36, No. 7 (1967). 

[27] D. V. Morgan and D. Van Vliet, Can. J. Phys. 46, 503 (1963); D. V. Morgan and D. Van Vliet, 
Radiat. Effects and Defects in Solids 8, 51 (1971). 

[28] D. van Vliet in "Channeling", ed. by D. V. Morgan (Wiley, London), 37 (1973). 

[29] J.U. Andersen and L.C. Feldman, Phys. Rev. B 1, 2063 (1970). 

[30] K. Komaki and F. Fujimoto, Phys. Stat. Sol. (a) 2, 875 (1970). 

[31] Chapter 2 of "Ion Implantation", by Geoffrey Dearnaley, Amsterdam, North-Holland Pub. Co.; 
New York, American Elsevier, 1973. 

[32] B. R. Appleton and G. Foti, "Channeling" in Ion Beam Handbook for Material Analysis, edited 
by J. W. Mayer and E. Rimini (Academic, New York), p. 67 (1977). 

[33] Chapter 4 of "Ion-induced electron emission from crystalline solids" , by Hiroshi Kudo, Springer 
Tracts in Modern Physics, 2002. 

[34] V. V. Rozhkov and S. V. Dyuldya, Pisma Zh. Tekh. Fiz. 10, 1181 (1984) [Sov. Tech. Phys. Lett. 

10, 499 (1984)]. 

[35] S. 1. Matyukhin, Teehnieal Physics 53, 1578 (2008). 

[36] K. Lonsdale, Acta Cryst. 1, 142 (1948); "X-Ray Diffraction" by B. E. Warren, Addison Wesley, 
p. 61 (1969). 

[37] A. V. Sharko and A. A. Botaki, Sov. Phys. Jour. 14, 330 (1971). 

[38] A. V. Sharko and A. A. Botaki, Sov. Phys. Jour. 14, 765 (1971). 

[39] N. Neelakanda Pillai, C.K. Mahadevan, Physica B 403, 2168 (2008). 

[40] P. Geeta Krishna, K. G. Subhadra and D. B. Sirdeshmuk, Acta Cryst. A54, 253 (1998). 

[41] A. W. Hewat, J. Phys. C: Sohd State Phys. 5, 1309 (1972); M. Inagaki, M. Toyoda, M. Sakai, 
Jour, of Materials Sci. 22, 3459 (1987). 

[42] K. M. Gorski et at, ApJ 622, 759 (2005). 

[43] J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961). 

[44] K. Komaki el al, Phys. Stat. Sol. (a) 4, 495 (1971). 

[45] N. Bozorgnia, G. Gelmini and P. Gondolo, "Channeling in direct dark matter detection IV: 
daily modulation of the WIMP signal" , in preparation. 

[46] N. Bozorgnia, G. Gelmini and P. Gondolo, arXiv:1008.3676 [astro-ph.CO]. 



-37- 



